1. Field of the Invention
The present invention relates generally to simulation of radiation transport. Particularly, the present invention relates to an improvement in radiation transport simulation. More specifically, the present invention relates to a radiation transport simulation system, program product, and related methods resulting in variance reduction.
2. Description of the Related Art
Regardless of which technique is used at the time of the diagnostic study to develop the radiation therapy treatment plan, in the delivery of either conformal radiation therapy treatments or static radiation therapy treatments, an accurate determination of radiation dose to the delivered is very important. Successful radiation therapy depends on accurately placing the proper amount of radiation upon the target without unnecessarily damaging surrounding tissue. Thus, it is necessary to relate the radiation dose determined to be delivered to the target at the time of the diagnostic study to how the radiation dose actually delivered at the time of the radiation therapy treatment. If the actual dose is not correct, the radiation dose may not be delivered to the correct location within the patient's body, possibly under-treating the target tumor or lesion and damaging healthy surrounding tissue and organs.
The Monte Carlo dose calculation method is generally considered the most accurate method to determine radiation or particle transport in heterogeneous media. Monte Carlo method has had many areas of application. For example, simulations of the physical processes of transporting neutrons and gamma rays through thick (meters) walls made of concrete and metals have been used in nuclear reactor designs, radiation protection and other purposes. The major obstacle of using Monte Carlo simulations is the slow computation speed resulting from the need to simulate tens of millions of particles at the region of interests. Depending upon various conditions, billions or even trillions of particles, depending upon the application, may need to be simulated from the radiation source. The majority of the simulation time is spent on the physical processes in the media between the radiation delivery apparatus and the region of interest. Most of the particles are absorbed and never reach the region of interest.
The efficiency of Monte Carlo code, ε, can be defined as:ε=1/(s2τ)where s2 is the variance and τ is the computing time. In the Monte Carlo simulation of radiation transport, each particle is called a history. If N particles are generated to represent the incident photon fluence, N histories in the medium will be followed. The values of s2 and τ are proportional to N−1 and N, respectively. For a given Monte Carlo code, ε is a constant and its numerical value depends on the software algorithm and the computer hardware.
With respect to radiation therapy, primary photon interactions in tissue, or phantoms as tissue-equivalent medium, effectively remove primary photons from the incident radiation beam, resulting in exponential attenuation of the primary photon fluence with the penetration depth z:Φ(z)=Φ0×exp(−μz)where Φ0 is the primary photon fluence at the phantom surface (z=0), Φ(z) is the sane at the depth of penetration z>0 and μ is the linear attenuation coefficient for primary radiation in the medium of interest.
With respect to neutron transport, the mechanisms of interaction are essentially different from photons. They can be classified into elastic interactions and inelastic interactions. During elastic interactions, neutrons elastically deliver energy to the nuclei (especially to protons) that make the medium. During inelastic interactions neutrons transfer energy to the nuclei, leaving them in excited states, which lead to the emission of photons. There are inelastic interactions such as nuclear reactions when neutrons are absorbed, which end up with the emission of photons, protons, etc. Finally, a non-interacting neutron decays into a proton and an electron (and an anti-neutrino) with a half-life of about 15 minutes. All the above mentioned interactions can lead to an exponential decay of the primary neutron fluence.Φ(z)=Φ0×exp(−μz)where Φ0 is the primary neutron fluence at the phantom surface (z=0), Φ(z) is the same at the depth of penetration z>0 and μ is the linear attenuation coefficient for primary radiation in the medium of interest.
For the Monte Carlo simulation process this means that the number of available particle histories N reduces with phantom depth z as:N(z)=N0×exp(−μz)
Because s2∝N−1, it can be seen that the variance of the simulated quantity, e.g. absorbed dose, increases as z increases:s2(z)=s2(0)×exp(μz)=k×N0−1×exp(μz),where k is a proportionality coefficient and N0 is the initial number of particle histories representing Φ0. To resolve this problem of deteriorating accuracy at large phantom depths z, i.e., to maintain the accuracy criterion, one can initially generate a larger number of particle histories N0, which inevitably results in longer computing time.
This is where variance reduction techniques become useful. Variance reduction techniques have been either mathematical approximations or “tricks” designed to reduce the computing time τ without increasing the variance s2 to thereby increase the efficiency of the Monte Carlo code. To evaluate a certain variance reduction technique, one can simply simulate the same case using Monte Carlo code without variance reduction and then with variance reduction implemented. When doing so, one requires (selects) a specific variance in the region of interest for both simulations. The efficiency gain of the variance reduction (VR) can then be defined as:Efficiency gain=CPU time (VR off)/CPU time (VR on).
The principle behind most of these techniques is that higher priority is given to fewer number of selected histories. These histories are the most important ones contributing to the dose in the region of interest. One strategy can include avoiding spending time to propagate particles contributing little to dose, or alternatively paying more attention to “useful” particles. This includes techniques such as, for example, Russian roulette and particle splitting, energy and range rejection, Kerma approximation, interaction forcing, and biasing toward under-sampled quantities such as certain scattering angles. Another strategy can include using analytical approximations whenever possible, especially during pre- and post-simulation processing. This includes techniques such as, for example, first collision, ray tracing, correlated sampling, and wavelet de-noising. A further strategy can include a particle re-using method. Such method calculates only a few samples using Monte Carlo, and then scales the results to different location or angles. Regardless of the strategies employed, it is recognized that the Monte Carlo results should be statistically unbiased. That is, the employed strategy should not distort the expected results.
Techniques that have been applied to photon transport in Monte Carlo simulations include, for example, interaction forcing, particle splitting, and Russian roulette. In interaction forcing, many photons pass through the medium without interacting with it. Time is spent tracking these photons even though they do not contribute to the dose. These photons can be forced to interact within the simulation geometry and contribute to the statistics. In particle splitting, when a photon approaches the region of interest, it is split into ns identical sub-photons, each carrying a weight factor 1/ns. The calculated dose is thus unbiased and less time is wasted on photons traveling outside the region of interest. In Russian roulette, when a photon moves away from the region of interest, it is “killed” with a certain probability, K<1, and, if it “survives”, its weight is increased by a factor 1/(1−K). Both particle splitting and Russian roulette leave the Monte Carlo simulation unbiased while computing time is spent more efficiently.
In a first collisions method, the time-consuming simulations of the primary photon interactions can be replaced with the analytically calculated collision density. This technique depends on the accuracy of pre-calculated attenuation rate of primary photons tissue. This technique is closely related to “ray tracing.”
U.S. Pat. No. 6,714,620, by Caflisch et al. titled “Radiation Therapy Treatment Method” describes variance reduction through particle splitting and range rejection. U.S. Pat. No. 6,301,329, by Surridge titled “Treatment Planning Method and Apparatus for Radiation Therapy” describes variance reduction techniques including a kernel approach characterized by a filter to reduce statistical noise, and methods of splitting particles in regions of interests while discarding most of them when exiting regions of interests. U.S. Pat. No. 6,772,136, Kant et al. titled “System and Method for Financial Instrument Modeling and Using Monte Carlo Simulation” describes variance reduction techniques through importance sampling used in financial modeling. U.S. Pat. No. 6,381,586, by Glasserman et al. titled “Pricing of Options Using Importance Sampling and Stratification/Quasi-Monte Carlo” also describes variance reduction techniques including importance sampling and stratified sampling. U.S. Pat. No. 6,518,579, by Xu et al. titled “Non-Destructive in-Situ Method and Apparatus for Determining Radionuclide Depth in Media” describes a variance reduction technique that biases photon emitting angles to maximize the usage of simulated photons. U.S. Patent Application No. 20030204126 by Rivard titled “Dosimetry for Californium-252 (252Cf) Neutron-Emitting Brachytherapy Sources and Encapsulation, Storage, and Clinical Delivery Thereof” also describes a method to bias particle emissions toward the region of interest.
U.S. Pat. No. 6,366,873, by Beardmore et al. titled “Dopant Profile Modeling by Rare Event Enhanced Domain-Following Molecular Dynamics” describes a method to enhance rare events when simulating ion implantations. In this variance reduction technique, an ion is split into two and their respective weights are decreased by half, hence the overall weight factors decrease with depth while the total number of ions remain constant with depth. U.S. Pat. No. 6,148,272, by Bergstrom et al. titled “System and Method for Radiation Dose Calculation within Sub-Volumes of a Monte Carlo Based Particle Transport Grid” describes a variance reduction technique that selectively tracks those particles most likely to pass through the region of interest.
A combination of variance reduction techniques may result in noticeable improvements in Monte Carlo efficiency. It has been reported that the combination of interaction forcing, particle splitting and Russian roulette increases the efficiency of Monte Carlo simulation by a factor of 3-4. The first collision technique is reported to yield a 2-fold improvement. There is no general “recipe,” however, as to which combination will perform the best.
Thus, recognized by the Applicants is the need for a new variance reduction system, program product, and related methods applicable to both homogeneous and heterogeneous absorbing media that require a significantly smaller number of initial photon histories to result in the same or better accuracy at the same selected depth of interest, that can maintain a primary particle (e.g., photon) fluence invariant with depth in the absorbing medium, and that can compensate for the artificial constancy of the particle fluence to yield unbiased results for simulated absorbed dose. Also recognized by the Applicants is the need for a new variance reduction technique which can be combined with other techniques to further improve simulation efficiency and/or accuracy.